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We present first-principles calculations of optimally lo- 
calized Wannier functions for Cu and use these for an ab- 
initio determination of Hubbard (Coulomb) matrix elements. 
We use a standard linearized muffin-tin orbital calculation in 
the atomic-sphere approximation (LMTO-ASA) to calculate 
Bloch functions, and from these determine maximally local- 
ized Wannier functions using a method proposed by Marzari 
and Vanderbilt. The resulting functions were highly localized, 
with greater than 89% of the norm of the function within the 
central site for the occupied Wannier states. Two methods for 
calculating Coulomb matrix elements from Wannier functions 
are presented and applied to fee Cu. For the unscreened on- 
site Hubbard U for the Cu 3d-bands we have obtained about 
25eV. These results are also compared with results obtained 
from a constrained local-density approximation (LDA) calcu- 
lation. 



I. INTRODUCTION 

During the past few decades powerful numerical 
methods have been developed for the ab-initio (first- 
principles) calculation of the electronic ground-state 
properties of solids. In most of these methods density- 
functional theoryEI (DFT) has been used to treat the 
electron-electron Coulomb repulsion, and a local-density 
approximation^] (LDA) (or a local spin-density approxi- 
mation, LSDA, for magnetic systems) has been used for 
the exchange-correlation potential. This procedure has 
been very successful for many materials and ground-state 
properties (e.g., crystal structure, lattice constant, bind- 
ing energy, and ionization energy), but it has its limita- 
tions; the band gap of semiconductors is not properly re- 
produced, for instance. Furthermore, for systems such as 
high-temperature superconductors, heavy fermion mate- 
rials, transition- metal oxides, and 3d itinerant magnets, 
i.e., for systems in which the Fermi level falls into a re- 
gion of narrow energy bands, the LDA is usually not 
sufficient. It is generally accepted that the problem for 
these materials is the strong electronic correlations that 
are responsible for their electronic properties. For a de- 
scription of such strongly correlated systems one usually 
instead uses as a starting point model Hamiltonians like 
the Hubbard modelj and its multi-band generalizations. 
But in these models the Coulomb (interaction) matrix 



elements and also the one-particle (hopping) matrix el- 
ements that determine the unperturbed band structure 
are usually treated as free, adjustable parameters, i.e., 
they are not known from "first principles" for the given 
material; on the other hand, Coulomb correlations can be 
studied within reliable many-body approximations that 
go beyond the Hartree-Fock approximation. 

Both the ab-initio LDA and the many-body model- 
Hamiltonian methods based on Hubbard-like models 
have their merits, but until rather recently they have 
been almost separate and complementary approaches. 
But, in view of the power of each, a combination of 
these methods is desirable; and, in fact, during the 
last few weass there have been some attempts in this 
direction.aii-3 All of these recent developments add lo- 
cal, screened Coulomb (Hubbard) correlations U between 
localized orbitals to the one-particle part of the Hamil- 
tonian obtained from an ab-initio LDA band-structure 
calculation, but differ in how they handle the correlcu 
tion part. In the earliest attempts, the LDA+TJ methodcl 
used essentially a static mean-field-like (or Hubbard-I- 
like) approximation for the correlation. The simplest 
approximation beyond Hartree-Fock, secaaiLarder per- 
turbation theory (SOPT) in U, was usedBOErlij to study 
the electronic properties of 3d-systems (like Fe and Ni) 
and heaw—fipmion systems (like UPts). The LDA++ 
approacroEjiid has a similar strategy, but uses other 
many-body approximations to treat the correlation prob- 
lem, namely, either the fluctuation exchange approxij 
mation (FLEX) or the dynamical mean field theory 
(DMFT). Some of the other many-body treatmentsQliilE-! 
have also used DMFT, which is based on the limit of 
large-dimension Xd — > oo) approximation for correlated 
lattice electrons.^ Within DMFT the selfenergy becomes 
local, i.e., independent of momentum k, which allows 
a mapping of the lattice problem onto an effective im- 
purity model. The LDA+DMFT treatmentsEm men- 
tioned above differ in the many-body method they used 
for the effective impurity problem, namely, Quantum 
Monte Carlco (QMC), the non-crossing approximatior£3 
(NCA) or (iterated) perturbation theoryQ (IPT). But all 
these approaches, including the LDA+U, have in com- 
mon that they have to introduce a Hubbard U as an ad- 
ditional parameter, and hence are not real first-principles 
(ab-initio) treatments. Although they use an LDA ab- 
initio method to obtain a realistic band structure, i.e., 
single-particle properties, Coulomb matrix elements for 
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any particular material are not known, and the Hubbard 
U remains an adjustable parameter. 

One can obtain estimates on the magnitude of U ei- 
ther from experiment (high-cnc^sy. spectroscopy) or from 
the constrained LDA method.EHlZrE3 Within the latter 
method one adds the constraint that the electron occu- 
pation number for the correlated bands is fixed to a given 
number through a Lagrange parameter. One can then 
use LDA to calculate the ground-state energy for differ- 
ent occupations of the correlated states, and the differ- 
ence between the energy for double and single occupation 
is an estimate for the Hubbard U. This method has the 
advantage that effects of screening are already somehow 
included. On the other hand there are usually several 
bands and many interaction matrix elements (on-site, 
density-density, intra-band, inter-band, exchange, inter- 
site, etc.) that have different magnitudes, and the con- 
strained LDA can only give some average value for these 
various Coulomb matrix elements and not the individual 
ab-initio parameters (Coulomb matrix elements). This 
approach is intuitive and contains some type of screen- 
ing within a one-electron LDA approach. It is difficult to 
sort out the actual approximation involved. 

In this paper we suggest a different approach, namely, 
the direct ab-initio calculation of the one-particle (tight- 
binding) and two-particle (Coulomb) matrix elements. 
Our starting point is a standard electronic band-structure 
calculation, for which we h ava u sed the linearized muffin- 
tin orbital (LMTO) methodt§n.a within the atomic sphere 
approximation (ASA). The LDA band-structure calcu- 
lation yields not only one-particle energies but also 
their eigenstates, the Bloch wavef unctions, which form 
a proper basis of a one-particle Hilbert space. To deter- 
mine the local (on- and inter-site) Coulomb matrix, e b e- . 
ments it is necessary to construct Wannier functions,c3'E£l 
which are closely related to the Bloch functions via a uni- 
tary transformation, but which are not unique since the 
phases of the Bloch functions are undetermined. _ 

As first suggested by Marzari and Vanderbilt,cZl this 
gauge freedom can be used to construct "maximally lo- 
calized Wannier functions." Those are just Wannier func- 
tions with a special gauge that makes them optimally 
localized according to some criterion. A proper localiza- 
tion of the Wannier functions is important in our opin- 
ion, because only then do the standard assumptions of 
the model treatments hold such that only a few (on-site, 
nearest, and next-nearest neighbor) one-particle (hop- 
ping) and two-particle (Coulomb) matrix elements have 
to be considered explicitely. These matrix elements can 
then be calculated from the Wannier functions. 

We use two different methods to calculate Coulomb 
matrix elements from Wannier functions. The first 
method uses the fact that the LMTO-method provides 
Bloch functions in the basis of linear muffin-tin orbitals.Eil 
Therefore the Wannier functions are given as linear com- 
binations of such muffin-tin orbitals as well, and can be 
used to evaluate the Coulomb integrals efficiently, simi- 
larly to what was done in Ref. E8[ The second method 



uses a fast Fourier transformation (FFT). It does not rely 
on the property of the wave functions being linear and is 
therefore more general. It is also very quick and efficient. 

The paper is organized as follows. In Sec. || we present 
some of the computational details: we describe the form 
of the Bloch functions in the LMTO method, how to 
obtain the Wannier functions from them, and how to op- 
timize the choice of the Wannier functions by using the 
Marzari- Vanderbilt method. Then we describe how the 
one-particle (hopping) matrix elements are obtained from 
these localized Wannier functions, and we present the 
two methods to calculate the Coulomb matrix elements. 
To illustrate the method we have performed actual cal- 
culations for a well understood system, namely for Cu. 
Although this material is not a strongly correlated sys- 
tem, it has almost completely filled, narrow, 3d-bands, 
for which (well localized) Wannier functions and one- and 
two-particle matrix elements can be calculated. Results 
for Cu are presented in Sec. [II, where we show some of 
the Wannier functions, demonstrate how well localized 
they are and that the one-particle (tight-binding) matrix 
elements obtained from them allow for a reconstruction 
of the band structure. The direct Coulomb matrix ele- 
ments obtained are rather large, between 20 and 25 eV for 
Wannier states with mainly 3d-character, and are about 
5 eV for nearest neighbors, and about 1 eV for exchange 
interactions. In Sec. |^we describe constrained LDA cal- 
culations, which yield somewhat smaller values (about 18 
eV) for the Hubbard U of Cu, and in the final section (|y|) 
we discuss how to extend and further apply the current 
approach. 



II. COMPUTATIONAL DETAILS 

We restrict ourself to the case where there is only one 
atom per unit cell and where we can neglect spin (non- 
spin-polarized calculations). For a given material the 
only input to an ab-initio LDA calculation i^-ithe atomic 
number. In the density-functional approach,EH a local- 
density approximation is normally used for the exchange 
and correlation interactions between the electrons; we 
have used the von Barth-HedintHI exchange-correlation 
potential and a frozen-core approximation. Within DFT 
the total energy of the ground state could be calculated 
as a function of volume for a given crystal structure and 
used to determine the equilibrium lattice constant; it is 
usually in good agreement with experiment. However, 
since our focus is on the determination of Coulomb ma- 
trix elements in a Wannier basis, we have simply used 
the experimental lattice parameters. 



A. LMTO wave functions 

For ouEJaand-structure results we have used the LMTO 
methocSEj within the atomic sphere approximation 
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(ASA). The combined correction termed was not in- 
cluded. The muffin-tin spheres are overlapping and their 
radius (the Wigner-Seitz radius S) is determined by the 
condition that the sphere volume equals the volume of 
the unit cell. Within the muffin-tin spheres the potential 
and wave functions are expanded in spherical harmon- 
ics with a cutoff Z max = 3, i.e., s, p, d, and f-orbitals 
are included. Furthermore the Bloch wave functions are 
given in terms of the solution to the radial Schrodinger 
equation (j> v i (r) to some fixed energies E v \ and its energy 
derivative <j> v i(r): 

*»k(r) = J2 {Mr)At + j>ui{r)Bt) Y L (r) , (f ) 

L 

where we use complex spherical harmonics in all of our 
calculations. This expansion is valid in one muffin-tin 
sphere. Here, as usual, L = {I, m} is understood and n is 
the band index and k is the wave vector. We define n by 
the condition that E n (k) < E n+ i(k). The virtue of using 
this method for Wannier functions is the simplification 
that only integrals over spheres are needed; no real-space 
integrations over complicated Wigner-Seitz (unit) cells 
are required. 

The Bloch functions obey 

* nk (r + R)=e ikH * nk (r) . (2) 

Therefore, the knowledge of a Bloch function in a single 
muffin-tin sphere is sufficient for the knowledge of the 
function in the whole crystal. This situation is differ- 
ent when we consider Wannier functions, which can be 
centered on different sites. It is useful to introduce a no- 
tation that holds for both Bloch and Wannier functions. 
To do this we perform an expansion like Eq. (|f|) in each 
muffin-tin sphere, which we label by its site vector R. 
The complete wave function (either Bloch or Wannier) is 
then given by: 

* a (r) = ^$ a (R i ;r-R i ) (3) 

i 

In this equation we have used the general notation for the 
wave function expansion $ a (Rj; r — R,-) such that: (i) <I> 
is any kind of wave function, (ii) The a stands for quan- 
tum numbers (Bloch: a — {n, k}; Wannier: a = {R, n}). 
(iii) The first argument in the parenthesis indicates the 
muffin-tin sphere about which we are expanding and is la- 
beled by its site vector, (iv) The second argument in the 
parenthesis is the position inside this muffin-tin sphere 
described by its relative vector. This means that this vec- 
tor has zero length in the center of the muffin-tin sphere 
described by the first argument, (v) Note that, for every 
R, 

i> a (R;r)=0 if|r|>S\ (4) 
In the case where $ is a Bloch function we find 

tf„ k (R;r) = e ikR *„ k (r) . (5) 



It is easy to see that Eq. (||) inserted in Eq. (||) obeys 
Eq. §. 

Also note that Eq. (|^) disregards the effects of over- 
lapping muffin-tin spheres. Within the ASA approxima- 
tion, all derivations are done as though non-overlapping 
muffin-tins are being used, and then these formulas are 
used with expanded muffin-tins, whose volumes sum to 
equal the unit cell volume (where the muffin-tin radius is 
expanded to a Wigner-Seitz radius for one atom per unit 
cell). In addition, this approximate eliminates the ne- 
cessity to handle interstitial regions, and hence the ASA 
formalism is mathematically much simpler than a full- 
potential electronic-structure calculation would require. 

B. Wannier functions 

In this section we show how to calculate Wannier func- 
tions from the LMTO tape of Bloch functions of Eq. (Q). 
The Wannier functionsE3 are defined by 

w Rn (r) = (r|Rn) = 1 £ e-* kR * nk (r) . (6) 

k 

Here, N is the number of k-mesh points in the Brillouin 
zone or, equivalently, the number of unit cells in the 
real space supercell that is used to discretize the k-mesh. 
As mentioned above, Wannier functions are not unique. 
Consider, for example, a single band n with Bloch func- 
tions l^nic); a transformation of the kind 

|*„k> - e i0 »'|*„k) , <4 k) real, (7) 

will still lead to Bloch functions. We shall call this a 
gauge transformation o£.the first kind. In the case of a 
composite set of bands JH3 this nonuniqueness corresponds 
to the freedom to choose the phases and "band-index 
labeling" at each k point of the Bloch functions: 

|*«k> - £te2l* m k) (8) 

m 

We shall call this a gauge transformation of the second 
kind. Here Umn is a unitary matrix. From all the arbi- 
trary choices of Wannier functions we will pick out that 
particular set that minimizes the total spread given by 

n = J2[(r 2 )n-(r)i] . (9) 

n 

For any operator A, (A) n denotes the expectation value 
(Rn|A|Rn). A method for minimizing-jEq. (||) has been 
developed by Marzari and VanderbilinJ and its applica- 
tion to the ASA wave functions does not pose any par- 
ticular problems (details will be given below). 

Before minimizing Q according to this procedure, it is 
useful to prepare the Bloch orbitals to make the start- 
ing Wannier functions somewhat localized. This has two 
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advantages: (i) the minimization procedure converges 
faster, and (ii) this helps avoid gettisg trapped in lo- 
cal minima. Marzari and VanderbiltE3 suggest several 
possible preparations. We have found our own method, 
which seems to work well. This involves a simple gauge 
transformation for each band, which is given by 



*nk(r) -> exp(-ilm ln*„ k (r ))* nk (r) 



(10) 



This gauge transformation has the property that 
Im In ^ , „k(ro) transforms to zero. So at the point ro 
all the Bloch functions will have the same phase (in this 
case just 1 + iO) and (ro|On) will take a large value. We 
thus expect the Wannier function to be fairly localized 
at ro. To make the method work well one should choose 
Tq where the Wannier functions are expected to be rea- 
sonably large. In our calculations we have chosen the di- 
rection of this vector to be well away from the expected 
zeroes of the spherical harmonics and with an absolute 
value far enough away from the origin to be in a place 
where the Wannier functions should have a significant 
magnitude. We found an ro of (0.8, 1.0, 0.3)ao to work 
well for fee Cu. 

We shall now derive expressions of the form of Eq. (||) 
for Wannier functions. From Eqs. (o) and (pi) we have 



( R '; r ) = ^E e " kRvI/ «k(R';r) 

k 

= ^E elk(R " R) *«k(r) 



(11) 



Because Wannier functions on different sites have the 
same form (shape) of their wave functions and differ only 
by a translation of their origin, it is useful to use a nota- 
tion that indicates values of a wave function relative to 
a Wannier function centered at the origin: 

w Rn (R'; r) = iz> „(R' - R; r) = iy n (R' - R; r) , (12) 

where we have introduced the notation wo n = w n (i.e., if 
the subscript contains only a wave function label without 
a spatial vector R, then we are using a relative notation 
that refers to a Wannier function centered at the origin) . 
We can use the Bloch condition (cf. Eq. (pd|)) to calculate 
the parts of the Wannier function on other sites R: 



w„(R;r) = jrJ2' 

k 



ikR 



*nk(r) . 



(13) 



Note that |r| < S. For this notation to work in our 
numerical calculations, it is essential to force the Wannier 
center, i.e., the muffin-tin sphere where (r|0n) is largest, 
to be at the muffin-tin sphere around the lattice site 0; 
we achieve this by setting r | < S in Eq. jl0|). In most 
of the rest of the paper, we will almost always use the 
relative notation that refers to Wannier states centered 
at the origin, and will perform whatever translations are 
necessary to be able to use these states. 



In the method of Marzari and VanderbiltJ^ the start- 
ing point for the calculations are a set of reference ma- 
trices defined by 



Mi°)( k > b > = (* mk | e -' ibr |*„. k+ b) • 



(14) 



Here b denotes a nearest-neighbor vector on the dis- 
cretized mesh in k-space (in this method, the set of b- 
vectors are needed for numerical derivatives). We calcu- 
lated the action of e~ jbr on the ket by using Eqs. ( A4) 
and (|A2[) and solved the remaining integral by using 
Eq. Q22[). We used a uniform (cubic) discrete k-mesh 
with a spacing Ak of 0.2(27r/a). In such a mesh there 
are 6 nearest-neighbors for the b vectors needed for the 
numerical derivatives. We were careful not to double 
count vectors in the k-mesh (those equivalent to each 
other by a reciprocal lattice vector) within the Brillouin 
zone (which has 500 k points in the full zone). 

We then used the steepest-descent method and rele- 
vant equations in Sec. IV of Ref. 27 to iterate a series 
of small steps where a set of AW^ were calculated and 
used to update the unitary matrices C/ k and the M k,b 
matrices. After each iteration, where we update all the 
relevant k matrices, we calculated the spread function 17, 
and continued iterating until this converged. 

In these calculations the initial matrices M(°)( k ' b ) are 
by far the most time consuming computationally (it re- 
quires storing 6 • 500 • 16 2 = 768, 000 complex numbers). 
The iterations of the steepest descent method were much 
faster. For this reason we used many iterations (about 
1500 steps) and converged fl to about 0.01%. For the 
step size (cf. Eq. (57) of Ref. ^) we used an a of 0.2. 

The final result can be written in a form similar to the 
LMTO wave functions, Eq. (0), 

Wn (R;r) = £ (Mr)Al R + ^/(r)S£ R ) Y L (r) , (15) 

L 

where the A and B matrices originally come from the 
LMTO wave functions, but are then updated from the 
relevant phase information, unitary matrix, and other 
integrations and transformations of the method. 

Because of the normalization of the starting LMTO 
Bloch wave functions (which are normalized to unity 
within a single unit cell) , each Wannier function is natu- 
rally normalized to unity when integrated over all space. 



C. One particle matrix elements 

The Wannier function basis can be viewed as an or- 
thogonal tight- binding basis. For this reason it is useful 
to calculate one-particle matrix elements of the Hamilto- 
nian in the Wannier basis. As we shall sec, these matrix 
elements are (for a gauge transformation of the first kind 
only) equivalent to the Fourier components of the band 
structure; this equivalence is useful for checking some of 
the numerical aspects of the calculations. 
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Because the Hamiltonian has the property that H (r) = 
H(r + R), it is sufficient to consider the matrix elements: 

t Rnrn ee (Rn\H\Om) (16) 

Inserting Eq. (|) and using H\^ nk ) = £;„(k)|*„ k ) we 
find 

tR nm = ^E e4kR ^( k )' ( 17 ) 
k 

which are just the Fourier components of the band struc- 
ture. The Bloch states |^ nk ) continue to be eigenstates 
of H under a gauge transformation of the first kind and 
one can also easily show that the tn nm are invariant un- 
der this type of gauge transformation. The tji nm from 
Eq. ( p"7j ) can be directly calculated from the band struc- 
ture E n (k). 

A gauge transformation of the second kind leads to 
states |\&nk) that are no longer eigenstates of H with 
eigenvalue E n (k). Therefore tji nm are not invariant un- 
der a gauge transformation of the second kind. However, 
generally we can always calculate 

H* m EE (* nk |ff|tf mk ) = e~ lkR tn.nm (18) 

R 

and use its diagonalized eigenvalues as the band structure 
for any set of Wannier functions. The matrix H^ m is 
Hermitian; it is, of course, already diagonal for a gauge 
transformation of the first kind, with the energy levels as 
the diagonal matrix elements. 

To calculate tn nm from Eq. Jl6|), i.e., using Wannier 
functions, we can use Eqs. (0) and ( |T^ ) to find 

t Rnm = J2f d3r w »( Rl - R; r ) Hw ™(^ r) , (19) 
i J 

where the integral is over a single sphere only. The effect 
of H on the second wave function can be carried out 
easily because we are working in a linear basis. We only 
note that (H — E v i)<f>i(r) = and {H — E v i)(f>i{r) = <£,(r), 
for details see Ref. |2^. In order to calculate Eq. (|l9|) , we 
must evaluate integrals of the form 

/ = J d 3 rf*(r)h(r) , (20) 

where the functions /i(r) are given by the expansion 

fi(r) =^ l R iL (r)Y L (T) . (21) 

L 

Inserting Eq. ( pl| ) into Eq. ( pp| ) and using the orthonor- 
mality of the spherical harmonics yields: 

/ = ]T / 'dr r 2 R* 1L (r) R 2L (r) (22) 

L J 



Because Eq. (gj) was our starting point, the radial func- 
tions R(r) will always be given in terms of <f>i{r) and 
<j>l(r), i.e., 

R lL {r) = AiLcPtir) + B lL Mr) . (23) 

We will use this form to calculate the integral I very 
efficiently. It is clear that any integral can be reduced 
to a linear combination of "basic" integrals. Those basic 
integrals consist of the (very limited) combinations of the 
4>i (r) 's and <pi (r) 's. We will label them by 

bi; PlP2 = J drr 2 [S Pl , (f>i(r) + 5 Pltl <f>t(r)} 

x[5 P2! o<f>i(r) + 5 P2t i<j>i(r)] , (24) 

where p\ and p 2 can take the values and 1. So it must 
be possible to write the integral / as: 

l l 

1 = X! X! a i;piP2 & ';piP2 (25) 

L pi=0p2=0 

It follows that the coefficients aL ;Pl p 2 are given by 

a L;p iP2 — L^pi,0^1i + <5pi,i-Bix] 

x[S P2 ,oA 2L + 5 P2A B 2L ] ■ (26) 

We are now in a position to calculate Eq. (^) with the 
aidofEq. @. 

D. Wannier-Function Projected Density of States 

The density of states (DOS) per spin is defined by 
AT(£) = _^W d 3 k 6(E - E n (k)) , (27) 

( Z7r ) ~ JBZ 

where V is the volume of the unit cell. In the same 
way that the DOS is often projected in terms of the l- 
character of the states, one can do a similar treatment 
for a projection onto the Wannier states. We can define 
a projected DOS for Wannier states, by inserting the 
projection operator onto the Wannier states |0j)(0j| into 
Eq. (H): 

Nj(E) = -^J2j BZ d " k \(^\0j)\ 2 S(E - E n (k)) 

(28) 

Note that the ijj n k in this formula have to be the Bloch 
states before the gauge transformation, since the band 
structure E n (k) is related to the untransformed states. 
The Bloch wave functions are normalized to a single unit 
cell and each Wannier function over all space. We, can 
calculate Nj(E) by using the tetrahedron method.Eil For 
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the k-points that form the tetrahedras we need to cal- 
culate |(V>nk|0.7'}| 2 , whic h we have done using the scheme 
described in section [II C| . In these calculations it is impor- 
tant to be aware that \0j) has parts of its wavefunction 
on other sites than the central site where it is centered. 
In our calculations, we included parts of the Wannicr 
function out to 17 near-neighbor shells of sites. 

Note that the exact projection operator is a sum over 
all R, since 



£|Rj>(Rj| 



1 



(29) 



However, it is sufficient to only consider the Wannicr 
states \0j) in our projection (and not all the |Rj)), since 



|<v„k|Rj')| = |(^„ k |oj)| 



(30) 



We can check the correctness of our projection by com- 
paring N tot (E) = £\ Nj(E) with the N(E) that is cal- 
culated directly from the LMTO energy eigenvalues. We 
find that our projected sum is accurate to within 0.2% of 
the LMTO value. 



E. Coulomb matrix elements 



The matrix elements we wish to calculate are 



W- 



12,34 



r — r 



—w^w^ws^w^r') . (31) 



where 1, 2, 3, 4 = i = {R^} is a Wannier state, and W 
denotes the Coulomb interaction. The spatial integrals 
over r and r' extend over all space. Using Eqs. (||) and 
(|l2|), we can use translations to rewrite this expression 
so that the integrals are only over the muffin-tin sphere 
at the origin: 



W 12M = Y,W{12,34;TL i ,R j ) 



(32) 



where the expression W(12, 34; R, R') is defined by: 

| r f r r /+ r R _ e2 R1 < ( R ~ R i; 'X ( R - R2 ' r ') 

x^„ 3 (R-R 3 ;r)iy n4 (R / -R4;r / ) (33) 

Since most applications of the Hubbard model use only 
two orbitals instead of all four, it is useful to define the 
limiting subset of the W functions as direct Coulomb fTy 
and exchange integrals: 



Ui2— Wi2,12 
Jl2= Vvi2,21 



(34) 



and the obvious generalizations for: 

[7(12; R, R')= W(12, 12; R, R) 

J (12; R, R') = W(12, 21; R, R) (35) 



1. Spherical harmonics expansion 

We will now only consider matrix elements between 
Wannier functions centered on the origin (i.e., where the 
Ri = in Eq. (|3l]). Because we are using maximally 
localized Wannier functions, most of the Wannier func- 
tions have their largest component in the center cell (see 
section III). As a first approximation, we will therefore 



neglect all other muffin-tin spheres. This approximation 
allows us to calculate on-site inter-band matrix elements. 
We are thus looking for 



W nin2m ~W (12, 34; 0,0) 

= J d 3 rdV<(0;r)<(0;r') 

e 2 

X | r _ r/ , Wn 3 (0; r)w ni (0; r ) , 



(36) 



where the integral over r is only over the central site. 
Inserting the expansion Eq. (21) for the Wannier func- 
tions and making use of the well-known expansion (see 
for example Ref. |3^) 



1=0 



where r > (r<) is the length of the greater (smaller) of 
the two vectors r and r', we find 

/ = E^TTE / drr*R* 1Ll (r)R 3L3 (r) f dr'r' 2 x 

l Li J 

I I 

x^ L2 (r')i?4L 4 (/)4fi E C ^ ■ (38) 



The coefficients Cll'L" are called Gaunt coefficients [see 
Eq. ([A3]) in the appendix] . If we define 



Cl;L 1 L 2 L 3 L 4 . — ^ Cl 3 L ± L C'l 2 L4.L 



(39) 



and 



drr'Rl Ll (r)R 3Ls (r) 



x / dr'r' Z R* 2L2 (r')R iLi (r')- rHl 



the integral takes the form 



E 

l, Li 



47T 

21 + 1 



G: 



(40) 



(41) 



The task is now to determine h-Li (we use the short- 
hand notation Li for LiL 2 L 3 L4). To do this, we will use 
the formalism developed in the last section. In complete 
analogy to Eqs. (p4|-|2q) we now find 



G 



where 



Pi 

2 
i=l 

4 
i=3 



(42) 



(43) 



and 



°l;hp t 



dr r 2 {5 Pl fi^ h {r) + S Pl ,i4>h(r)] 
x / dr' r' 2 [8 P2fi (f) l2 {r') + 5 P2A <j) h (r')} 



x[5 p4 , ^ 4 (r') + <5 P4) i4(^)]^Tr • ( 44 ) 

r > 

It should be noted that these basic integrals are symmet- 
ric with respect to some of their indices. If we introduce 
the joined index rii — {U,pi} then: 

^l^nin^n^n^ — ^l\n^U2Ti\n^ — ^;m7l4TO3n2 — 
^Z;n3n4ni7i2 — ^l\n^n\ri^n^, — ^l\nin^,ri^n\ — 
^i;n4nin2n3 — ^Z;n4n37i2ni (^^) 

If we consider the numerical aspects for the case where 
s, p, d, and / orbitals are included in the wave func- 
tion expansion, we find that we need to use a cutoff of 
'max = 6 in Eqs. ([37]) and (|39|). Using the symmetries in 
Eq. (|45|), we then find 9072 basic integrals that have to 
be calculated and stored. The sum in Eq. ([H]) however 
involves 7 ■ 16 4 = 458, 752 elements. Fortunately, only 
6778 combinations of the I, L\, L2, £3, L4 coefficients in 
Eq. (IH) have to be calculated; the others vanish. Each 
of these coefficients involves a sum over 16 elements, and 
each of these elements is a product of 5 numbers. 



2. Fast Fourier Transformation (FFT) approach 

The method we have just described works well, but 
requires a lot of Gaunt functions and other complications. 
As written, it also only involves integrals over the central 
site and ignores parts of the Wannier functions on nearby 
neighbors. We have therefore found a different approach 
to the problem. 

To calculate W(12, 34; R, R') for any lattice sites R 
and R', we make use of the Fourier transform 



— 5- = TT 



(46) 



PF(12,34;R,R) = ^J^ e ^ R ')/i 3 (q)/ 24 (-q) 
/y(q) = j d 3 re 4qr <(R " r)w nj (R - R,;r) (47) 

The fij functions are just the Fourier transforms of a 
product of some Wannier functions in a sphere. These 
can be calculated very efficiently by calculating the Wan- 
nier functions on a cubic mesh in real space and then 
applying a standard FFT algorithm. To do this, we have 
used the routine "fourn" (cf. Rcf. |33|). For details on 
how to apply the FFT to continuous functions, Ref. [54] is 
very useful. The result of the Fourier transform is /«(q) 
on a cubic mesh in q-space with some Aq (the distance 
between the mesh points). We perform the remaining 
q-integral in the following way. Let us call the integrand 
without the q~ 2 term 



F(q) = e^ R - R '>A 3 (q)/ 2 4(-q) 



(48) 



which is smooth function at q = 0. In order to treat the 
divergence arising from q~ 2 , we split the integral in the 
following way: 



F(0) 



d 3 q 

„2 



(49) 



All integrals are over a cube with length NAq. The last 
integral is just half of this length times C, which we define 

as 



C 



dx 



dy 



dz — 



15.34825 



(50) 



The remaining integral in Eq. ( [49] ) is transformed into a 
sum over little cubes with volume (Ag) 3 . For q = the 
value of integrand is calculated via the second derivative 
of -F(q) at q = numerically (the second derivative is 
needed to cancel the q 2 in a power-law expansion of F). 
The cubic grid in real space that we used to calculate 



the Wannier functions in Eq. 



had N = 64 3 points 



in the real space grid with a spacing Ax = 0.17. The Aq 
spacing of the g-mesh is determined by N and Ax. Using 
the FFT for continuous Fourier transformations one has 
to be very careful about the choice of these values because 
the FFT is a discrete Fourier transform. It is important 
to make sure that the results of a FFT calculation do not 
depend on the values N and Ax. 



and find for Eq. (B3| 



Note that each integral in Eq. (47) could be calcu- 
lated from the spherical-harmonic expansions. However, 
many such integrals would be required and the method 
would be extremely computationally expensive. The 
FFT method generates all the q values needed with a 
single calculation and is much more efficient. However, 
because of finite mesh sizes and compromises between 
real and q space integrals, it is not as acc urate as the 
spherical-harmonic expansion method of Sec. [I E 1 , when 
the latter is applicable. 
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III. RESULTS AND DISCUSSION 

We have tested our methods on Cu, which has the fol- 
lowing properties: (i) It has a simple close-packed foe 
crystal structure for which the ASA should be a rea- 
sonable approximation, (ii) Cu is a simple metal that 
belongs to the 3d transition metals, so one can deter- 
mine Coulomb matrix elements for the 3c? states which 
are interesting and of relevance for the really correlated 
3d-systems. (hi) Since Cu is non-magnetic, we do not 
have to worry about spin-polarized or magnetic calcula- 
tions. We have used the experimental lattice constant 
614A as given in Ref. [55]. As usual we use atomic 
Rydberg units and ao = h 2 /me 2 is the Bohr radius. 



5 10 15 |R| • a„ 




FIG. 1. Localization of Wannier functions before (grey) 
and after (solid) minimization of the spread functional f2. 
Each dot represents the portion of the wave function in a 
muffin-tin sphere. The best exponential fit to the decay is 
roughly e ,r with the values of 7 given within the figure. In 
Eq. ^ we set r = (0.8, 1.0, 0.3)a . 

From the one-electron Bloch wavefunctions, the Wan- 
nier functions were obtained using 500 k-points in the 
BZ. Since we have used a cutoff of Z max = 3, the LMTO 
method generates 16 bands. In the minimization pro- 
cedure all 16 bands were treated as a composite set of 
bands. To demonstrate the localization of the Wannier 
functions obtained, we have calculated 

(w n \w n ) R = J d 3 r|w„(R;r)| 2 , (51) 

which is the relative weight the Wannier function local- 
ized at the site has in the muffin-tin sphere centered 
around R. In Fig. [l] we have plotted for n = 2 the func- 
tion hi(w2\w2)n. as a function of |R|. Although the con- 
tribution to W2 appears to decrease exponentially with 
increasing distance from the central sphere when plot- 
ted in this way, i.e., our Wannier functions are expo- 
nentially localized, we actually get just as good a ht 
through the scatter of the data with a power-law de- 
pendence with a power of about -7. It is not easy to 
numerically decide whether the decay is an exponential 
or a power-law dependence, since our Wannier functions 
are ultimately periodic in a supercell determined by the 



Ak spacing of the discrete k mesh used to construct 
them. In either case, the Wannier function is highly lo- 
calized. The grey dots and line in the figure show the 
result when only the phase has been adjusted accord- 
ing to Eq. (|l0|); then the Wannier functions have a rela- 
tively smaller decay constant of 7 = O^a^ 1 . The black 
dots and line show the result after the full localization 
(minimization) procedure of Ref. ^ has been applied 
by minimizing the full set of all 16 bands considered; 
clearly a much better localization with a larger decay 
constant 7 = 1.14oq 1 has been achieved. When we tried 
to minimize a smaller subset of bands (5 bands instead 
of the full 16) the decay factor was in between the other 
two values, with 7 w .Tla^ 1 (not shown in the hgure). 



TABLE I. Angular character of Wannier functions in the 
center muffin-tin sphere (top) and first shell, i.e., 12 nearest 
neighbors (bottom). 



1 


n = 


1 


2 


3 


4 


5 


6 


(4)8 


.102 


.013 


.000 


.008 


.026 


.134 


.187 


(4)P 


.297 


.131 


.058 


.042 


.151 


.373 


.392 


(3)d 


.407 


.663 


.895 


.886 


.716 


.323 


.256 


(4)f 


.096 


.140 


.024 


.039 


.064 


.069 


.051 


E 


.902 


.947 


.977 


.975 


.956 


.899 


.886 


(4)8 


.009 


.004 


.002 


.002 


.004 


.009 


.011 


(4)P 


.019 


.008 


.004 


.004 


.008 


.020 


.023 


(3)d 


.052 


.029 


.012 


.013 


.024 


.053 


.060 


(4)f 


.012 


.008 


.003 


.003 


.006 


.012 


.013 


E 


.092 


.049 


.021 


.023 


.041 


.094 


.107 



Here we should note that the Wannier functions are not 
pure in terms of their Z-character. Table | shows the an- 
gular character in the center muffin-tin (MT) and the 
hrst shell for the first 7 Wannier states. We see that for 
the states with n—0 to n=A the d-character is largest 
which suggests to call these states d-like states yielding 5 
d-states per spin direction as expected. But among these 
states the d-character is highest (nearly 90 %) for the 
states n = 2 and 3. Table | also tells us how much of the 
state is found in the center muffin-tin. We see that the 
state n = 2 has 97.7% in the center MT and only 2.1% 
in the next shell demonstrating how well localized this 
Wannier function is. The Wannier functions correspond- 
ing to n = 0, n = 5, and n = 6 have considerable 4s- and 
4p-character, and n — 5 and n = 6 have the smallest 3d- 
character. But they are very well localized as well, having 
at least 88% of their total weight already within the cen- 
tral muffin-tin sphere. On the other hand, all Wannier 
states are mixed with respect to their Z-character, since 
the minimization procedure mixes all the I characters. 

Figure H shows a few radial averaged Wannier functions 
in their center MT. We should note that the peak of the 
states n = and 6 for r — > does not contribute very 
much to matrix elements because of the r 2 in Eq. ^2%). 
Figure ^| may also be qualitatively compared with the 
Wannier function of Cu in Ref. 
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i i i 


5 = 2.669a 



0.5 1.0 1.5 2.0 2.5 r ' a o 

FIG. 2. Radial averaged Wannier functions in the center 
muffin-tin sphere for various indices n. 



From these Wannier functions we have calculated the 
hopping matrix elements t-R, nm according to Eq. (19) 
These can be inserted into Eq 



( |18| ) in order to de- 
termine an effective orthogonalizcd (diagonal S-matrix) 
tight-binding representation. The matrices i? k are not 
diagonal because the unitary matrix which was used in 
the minimization of the Wannier functions scrambled the 
different bands. However, we can still diagonalize £f k 
for each fc-point and compare the eigenvalues with the 
original LDA band structure. The results are shown in 
Fig. ||, where we have cutoff the R-sum in Eq. (jl^) to 
include only and the 3 nearest shells, i.e., 43 sites. We 
have found that the decay of tn nm as a function of |R| 
is a lot faster than that of the Fourier components of the 
band structure, Eq. (|r?|). If we just take Eq. (O) and 
recalculate the band structure according to Eq. (181), the 
agreement is a lot worse (for the same number of sites 
in tji nn ). This can be understood in the following way: 
Labeling the bands according to E n (k) < E n+ i(k) is not 
"natural" , therefore at points in fc-space where two bands 
cross each other E n {k) has a kink. Those kinks have non- 
negligible Fourier components with large |R|, which our 
cutoff sets to zero. 



W -4 

-6 



r xwl r kxu 

FIG. 3. Comparison of LDA band structure (dashed line) 
and the diagonalized eigenvalues of Eq. ([l8|), where 3 shells 
in the lattice sum were included. The bands are relative to 
the Fermi level at OeV. 





FIG. 4. Comparison of LDA band structure (dashed line) 
and the diagonalized eigenvalues of Eq. (^), where 8 shells 
in the lattice sum were included. 

f| we have included 8 shells in the lattice sum 
i.e., 141 sites. As we can see the two curves 
agree even better. These calculations are a test of the 
quality of the Wannier functions, i.e., how well the ma- 
trix elements obtained from these Wannier functions re- 
produce the known band structure. 

With respect to the magnitude of the hopping matrix 
elements, the t 0nm are largest and provide information 
about the positions of the bands. For a next neighbor 
R the |iftnm| are of the order of 0.3eV for d-states (and 
leV for the state with n = 0). For larger R the hopping 
matrix elements are less than 0.15eV for d-states (and 
less than 0.5eV for the state with n = 0). 




-8 -7 -6 -5 -4 -3 -2-10 1 E/eV 
FIG. 5. Total electronic density of states, relative to the 
Fermi energy. We have used a Gaussian broadening of 
2 mRyd to remove the spikes inherent in the linear tetra- 
hedron method. 

Next, consider the projected DOS (PDOS). In—these 
calculations we have used the tetrahedron methocO (see 
also Ref. ^3|) with 200 k-points and 691 tetrahedras in 
the irreducible part of the Brillouin zone. In Fig. || we 
have plotted the total DOS which can also be found in 
the literature (see Ref. |37|). 
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-6 -5 -4 -3 -2 -1 S/eV 

FIG. 6. Projected DOS for Wannier states through 3. 




-6 -5 -4 -3 -2 -1 E/eV 

FIG. 7. Projected DOS for Wannier states 4 through 7. 

Figur es p] and [7] show the projected DOS according 
to Eq. (J2q) for the Wannier functions through 7. It 
is interesting that states j — 2 and 3 have very similar t 
character (cf. Table ||) but a very different projected DOS 
in Fig. ||, i.e., they are peaked at different energies, and 
emphasize different parts of the d band. Table O shows 
the projected density of states (actually the percentage 
of the DOS in different states) and the projected num- 
ber of states evaluated at the Fermi level, where the jth 
projected number of states is defined as 



nj(E F ) 



dE Nj (E) 



(52) 



This is just the number of electrons in the jth 
state. Every state could maximally be occu- 
pied with 2 electrons (one for each spin direction). 

TABLE II. Projected number of states at the Fermi en- 
ergy, i.e., nj{Ep). For j > 8 (the numbers not given) 
nj(E F ) < 0.25. The second line shows the percentage of 
the DOS at Fermi energy, i.e., 100 ■ N j {E F )/N{E F ). 



3 





1 


2 


3 


4 


5 


6 


7 


8 


nj{E F ) 


1.17 


1.31 


1.59 


1.80 


1.34 


.66 


.90 


.51 


.78 


% DOS 


4.7 


2.1 


11.1 


7.1 


15.4 


13.0 


4.6 


3.2 


27.3 



We next consider a calculation of the direct Coulomb 
integral Uu for a d-like orbital with itself. As discussed 
above, the Wannier function for n — 2 has nearly perfect 
d-character and to a good approximation we can consider 



only its contribution in the central MT, i.e., at site 0. 
What we then calculate is the on-site Coulomb matrix 
element between two electrons (with different spin be- 
cause of the Pauli principle) at the same site in the same 
Wannier state, i.e., essentially the Hubbard- U in it s orig- 
inal sense.cl The method described in Sec. II E 1 yields 
U(dd; 0, 0) =25.26 eV while the method from Sec. |IIE2 
yields 25.16 eV for this quantity. But with the second 
method we are able to calculate all the elements involv- 
ing tails of the Wannier functions in other muffin-tins 
in the double sum in Eq. (|32|). When we do this and 



include sites where Ri and Rj are nearest neighbors, 
we get Udd =25.51eV, which shows that including the 
portions of the Wannier function on neighboring sites 
is a rather small correcti on on U for such a strongly 
localized function. Table III show these quantities for 
the Wannier functions n — through 6. Going be- 
yond nearest neighbors would have even a smaller ef- 
fect. Therefore, one can truncate the sums over higher 
neighbor shells for the Coulomb matrix elements, which 
converge faster than for the hopping matrix elements. 
The reason for this is that the Coulomb integral in- 
volves a product of four wave functions, whereas the hop- 
ping matrix elements involve only two wave functions. 

TABLE III. Onsite FFT [/'s. In the first line (onsite-[7) we 
have only included Ri = Ri = in Eq. (^), i.e., U(jj; 0, 0). 
The second line (NN-C/) shows the same quantity, where we 
have included nearest neighbors for Ri and Rj . 



j 





1 


2 


3 


4 


5 


6 


onsite- U 
NN-J7 


14.82 
16.29 


19.05 
19.81 


25.16 
25.51 


25.49 
25.86 


20.79 
21.44 


13.81 
15.32 


13.22 
14.97 



The FFT approach allows us to calculate Coulomb ma- 
trix elements for Wannier functions centered on different 
sites. We have done this for the states a — {0, 2} and 
b — {R, 2} where R is a nearest neighbor of 0. Both are 
d-like states. In Eq. ( |32] ) we have again included nearest 
neighbors for R^ and Rj. The result is U a b=5.87 eV. 
The largest contribution in the sum is U(ab; 0, R) =5.66 
eV, which is the contribution arising from the two center 
spheres of states a and b. 

Our first method is most useful for calculating inter- 
band (on-site) Coulomb matrix elements when the states 
are so well localized that we can neglect the contribution 
from neighboring spheres. We have calculated both the 
direct Coulomb matrix elements U nm and the exchange 
integrals J nm for all n and m. Here the n, m just in- 
dicates the band and all Wannier states are centered at 
site 0. The results are given in Tables and for 
the first 7 bands. We see that the on-site intra-band 
Coulomb matrix elements are largest for the Wannier 
states n—2 and 3, which have almost pure d-character. 
We also note that all the direct Coulomb matrix elements 
U nm are rather large, while the exchange Coulomb ma- 
trix elements J nm with n ^ m are rather small (note 
that the diagonal terms for both Uu and Ju are identical 
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34|)). When we compare the diag- 
Jnn) from Tables IV and [v| with 



by definition, cf. Eq. 
onal elements U nn {= 

the first line of Table |II| we note relatively large differ- 
ences for the states ti=0, 5, and 6, which have large s 
and p character, as can be seen from Table || (e.g. 0.52eV 
for rt=0). This leads to a peak in their charge density 
near r = as we can see from Fig. || for n=0 and 6. 
For those states our FFT calculations had a numerical 
problem because our real-space grid was too large (with 
Ax = 0.17 do). But for 77=2, we do not have a peak 
at r = 0, and the FFT approach does an excellent job. 

TABLE IV. Inter-band on-site matrix elements: Unm in eV. 





771 = 


1 


2 


3 


4 


5 


6 


71 = 


14.30 


15.86 


17.86 


17.45 


16.31 


13.24 


12.77 


1 


15.86 


19.02 


21.36 


20.98 


19.33 


10. ZO 


1/1/10 
14. 4Z 


2 


17.86 


21.36 


25.26 


24.13 


22.30 


17.09 


15.86 


3 


17.45 


20.98 


24.13 


25.26 


21.88 


16.69 


16.12 


4 


16.31 


19.33 


22.30 


21.88 


20.70 


15.76 


14.78 


5 


13.24 


15.25 


17.09 


16.69 


15.76 


13.23 


12.32 


6 


12.77 


14.42 


15.86 


16.12 


14.78 


12.32 


12.43 


TABLE V. Intcr-band 


on-site matrix elements 


■ J nm 


in eV. 




771 = 


1 


2 


3 


4 


5 


6 


n = 


14.30 


0.91 


0.73 


0.76 


0.92 


0.98 


1.34 


1 


0.91 


19.02 


1.22 


0.84 


0.91 


1.43 


0.69 


2 


0.73 


1.22 


25.26 


0.92 


1.14 


0.95 


0.58 


3 


0.76 


0.84 


0.92 


25.26 


0.99 


0.71 


0.62 


4 


0.92 


0.91 


1.14 


0.99 


20.70 


1.20 


0.79 


5 


0.98 


1.43 


0.95 


0.71 


1.20 


13.23 


1.22 


6 


1.34 


0.69 


0.58 


0.62 


0.79 


1.22 


12.43 



The Hubbard-?/ clearly depends on the specific shape 
of the Wannier functions. Intuitively, one expects bigger 
C/'s for more localized orbitals. As an example of this, we 
have calculated a less localized Wannier state (by doing 
fewer steps in the minimization procedure) . In this case, 
the highest d-character state, which is almost pure d-like, 
has only 58.5% of its charge density in the center muffin- 
tin, and 95.4% within the first 3 shells. For those three 
shells we have used the FFT method to calculate all 43 2 
terms (|33"|). We find a U = 13.8eV for this less localized 
d-state. 

We should also note that most model calculations as- 
sume very localized, pure (in Z-character) Wannier func- 
tions. In particular, they often assume that LDA or 
some one-electron-like treatment is adequate for non-d 
and non-/ electron states, and that the only explicit cor- 
relations that need to be included are related to onsite (or 
sometimes also nearest-neighbor) Coulomb C/'s for the d 
(or /) states. It is also often implicitly assumed that the 
non-d and non-/ states have some screening contribution 
to the effective C/'s in the model Hamiltonian. These 
types of assumptions raise some difficulties for us to con- 
nect our treatment to the model Hamiltonians, since the 



orthogonalization properties and mixing necessary for lo- 
calizing our Wannier functions scramble the ^-character 
of the resulting orbitals. Hence, our effective C/'s do not 
have a pure d or / character (or s or p). Also, since 
we calculate C/'s for all of the orbitals, we are implicitly 
including correlation effects for all orbital (s and p as 
well as d and /), and however the C/'s in our treatment 
are ultimately screened in some many-body treatment, 
this screening may be different from that assumed in the 
model Hamiltonians. We may ultimately be forced to 
use some kind of projection of our orbitals onto pure 
^-character states in order to make appropriate identifi- 
cation between our types of states and more conventional 
model Hamiltonians. 



IV. CONSTRAINED LDA 

Finally, we have done a constrained LDA calculationEll 
to obtain an estimate for the Hubbard U . In this method 
the Hubbard U is defined as the Coulomb energy cost to 
place two (in our case d) electrons at the same site. This 
is 



U = E(N d + 1) + E(N d - 1) - 2E(N d ) 



(53) 



Here E(N d ) is the ground-state energy with N d d- 
electrons. If we consider this energy as a continuous 
function of Nd, where we constrain the value of Nd to 
be away from its minimized value, then the Hubbard U 
is given by: 



U dd = 



S 2 E(N d ) 



(54) 



This constraint, which fixes the total number of d- 
electrons to be Nd, can be taken into account by adding 
a Lagrange parameter v d to the total energy; i.e., the 
energy of the constrained system is given by 

E{N d ) = rnin[£{n(r)} + v d { J d 3 rn d {r) - N d }] . (55) 

Here E{n(r)} is the usual band-structure energy and 
n d (r) is the <i-electron density. On minimization the 
extra term in Eq. J55| ) leads to an additional constant 
potential, v d , in the Kohn-Sham equations, which acts 
only on the I = 2 angular momentum components of 
the wave function. Within the LMTO method, this is 
accomplished by adding a constant potential, v d , when 
solving the radial Schrodinger equation for I = 2, and 
then calculating the total energy as a function of v d . 
Since each value of v d changes the d occupation number, 
the final result can be written as E(N d ). This depen- 
dence is shown in Fig. || and can be accurately fitted by 
a parabola, SE = jU dd Nj, with U dd = 18.2 eV. This is 
of the same magnitude as the result obtained from the 
direct calculation of the Coulomb matrix elements, even 
though one might expect a smaller value because of the 
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screening effects that are believed to be included in this 
calculation. In our calculations we have only used a one- 
atom unit cell. If a larger unit cell is chosen, one could 
do a variety of additional constraints (e.g., changing the 
(^-occupation separately on two different atoms). Such 
calculations could attempt to sort out more details of ef- 
fective Hamiltonians (perhaps even two-particle parame- 
ters). However, such calculations would take our work in 
a different direction from what we are interested. Also, 
given the intuitive nature of the constrained method and 
the difficulties in fitting such a large parameter space, it 
is not clear how useful the resulting parameters would be 
or their uniqueness. 



U d =18.2eV 




Sn d 

-0.2 -0.1 0.1 0.2 

FIG. 8. The total energy as a function of the effective 
change in d charge. The line is a quadratic fit. 



V. CONCLUSIONS 

We have shown in this paper that ab-initio band- 
structure methods can be used for a first-principles cal- 
culation of well localized Wannier functions, which is 
achieved by— jusing a method proposed by Marzari and 
Vanderbilt.EZI From these localized Wannier functions the 
on-site and inter-site one-particle matrix elements of the 
Hamiltonian can be calculated. A good localization of 
the Wannier functions is needed to keep tight-binding 
(hopping) matrix elements restricted to a small num- 
ber of near neighbors. The Coulomb matrix elements 
within these localized Wannier states can also be cal- 
culated and are similarly only important between on- 
site and nearest-neighbor Wannier functions. The result 
is thus an electronic multi band Hamiltonian in second 
quantization with first-principles one- and two-particle 
matrix elements. The Hamiltonian is of the form of an 
extended multi-band Hubbard model but without ad- 
justable parameters; the parameters are directly calcu- 
lated for a given material. The only approximations still 
involved are the ones inherent to the ab-initio band struc- 
ture method used (e.g., the muffin-tin assumption, the 



ASA approximation, the choice of linearized orbitals in 
the LMTO, and the "frozen-core" approximation), and 
the truncation in the number of bands (states) per site 
that is explicitly taken into account (truncation of the 
Z-sum). The resulting multi-band Hamiltonian that in- 
cludes the Hubbard-[/ terms, of course, still has to be 
studied within a reliable many-body method or approx- 
imation, e.g, a multi-band version of the DMFT as in 
Refs. 0, [fij g and 

Our Cu calculations yield on-site direct Coulomb ma- 
trix elements ("Hubbard-U's") of the magnitude of 20 
eV for 3d-Wannier-states and inter-site (Hubbard-U's be- 
tween nearest-neighbors) values-pf 5 eV. This is the mag- 
nitude discussed already earliera and similar to those for 
atomic 3d-states. These U-values are much larger than 
commonly expected or used in model studies. Although 
our calculated Coulomb matrix elements are unscreened, 
the constrained LDA, which includes some screening ef- 
fects, gives comparable magnitudes for U. Since dynamic 
screening due to the mobile electrons in the outer shells 
(bands) is taken into account automatically when using 
an appropriate many-body method, e.g., a generalized 
random phase approximation (RPA), the only screening 
that should be included in a better theory is a static, 
short range screening by the inner core electrons. How- 
ever, the (atomic like) electronic states representing the 
inner ("frozen") core are well known, and it should be 
possible to calculate their screening contribution from a 
(generalized) static Lindhard theory. In future work we 
plan to examine the static screening due to the inner core 
states, an application of appropriate (multi-band) many- 
body methods, and the application to more strongly cor- 
related 3-d materials such as iron, cobalt, and nickel. 
Any treatment of screening will, of course, have to be 
very careful that screening effects are not double counted 
(once in the explicit screening and then a second time 
when the many-body Hamiltonian is solved). Finally, 
although any localized orbitals could be used as the ba- 
sis for a many-body treatment, the approach we have 
used (of constructing localized orbitals from LDA band 
states) has the advantage that these orbitals are a good 
basis set for any states without strong electron-electron 
correlations, since LDA is believed to be an accurate ap- 
proximation in this limit. We can hope that an addi- 
tional more explicit treatment of the strong correlations 
by a many-body theory will correct and improve on the 
LDA starting point. 
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APPENDIX A: SPHERICAL HARMONICS 
EXPANSIONS 

Any function A(r) within a (muffin-tin) sphere may be 
expanded in terms of spherical harmonics: 



A{v)=Y J A L {r)y L {i) 



(Al) 



If two functions A(r) and B(r) are given via their coef- 
ficients Ar,(r) and B[ J (r) 1 then the corresponding coeffi- 
cients Fr,{r) of the function Fir) — A{y)B(y) are given 
by: 



F L (r)= J2 A Ll {r)B L2 {r)C LlLL2 

Li,L/2 

The Gaunt coefficients Cll'L" are defined by 

c LL , L „ = J d 2 n Y L (n) r L *,(Q) y l ,,{q.) 

— $m" ,m' —m 



(A2) 



21' 



1 



-in 



J {L',L) 



(A3) 



and the c k (L',L) are tabulated in Ref. [58]. We may use 
Eq. (|A2|) to multiply a function with a plane wave e k (r) = 



-ikr 



whose coefficients are given by (see Ref. 



el(r) = AnMkr) \i l Y L (k) 



(A4) 
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